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The accuracy of two grid adaptation strategics, grid redistribution and 
local grid refinement, is examined by solving the 2-D Euler equations for the 
supersonic steady flow around a cylinder. Second- and fourth-order linear 
finite difference shock-capturing schemes, based on the Lax-Fricdrichs flux 
splitting, arc used to discretize the governing equations. The grid refine- 
ment study shows that for the second-order scheme, neither grid adaptation 
strategy improves the numerical solution accuracy compared to that cal- 
culated on a uniform grid with the same number of grid points. For the 
fourth-order scheme, the dominant first-order error component is reduced 
by the grid adaptation, while the design-order error component drasti- 
cally increases because of the grid nonuniformity. As a result, both grid 
adaptation techniques improve the numerical solution accuracy only on the 
coarsest mesh or on very fine grids that arc seldom found in practical appli- 
cations because of the computational cost involved. Similar error behavior 
has been obtained for the pressure integral across the shock. A simple 
analysis shows that both grid adaptation strategics arc not without penal- 
ties in the numerical solution accuracy. Based on these results, a new grid 
adaptation criterion for captured shocks is proposed. 


Key Words: high-order accuracy, grid adaptation, local grid refinement, grid redistribu- 
tion, shock capturing, finite difference scheme. 


1. INTRODUCTION 

Wave propagation phenomena in computational fluid dynamics, computational 
aeroacoust.ics, computational electromagnetics, and large eddy or direct numerical 
simulation of turbulence are characterized by the presence of both a lot of struc- 
ture in the smooth part of the solution and strong discontinuities. For this class of 
problems, attaining the design accuracy of high-order shock-capturing methods is 
problematic. As has recently been shown for 1-D unsteady and 2-D steady shocked 
flows, the first-order error introduced by shock capturing can persist globally down- 
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stream [1, 2]. The result is that the numerical solution is just first-order accurate 
downstream of the discontinuity, regardless of the design accuracy of the discretiza- 
tion used. Similar degeneration in accuracy for captured discontinuities has also 
been reported in [3]. 

One way of removing the first-order error component from the numerical solu- 
tion is shock fitting. As shown in [1, 2], if a cell interface is aligned with the shock, 
and a high-order conservative essentially nonoscillatory (ENO) formulation based 
on a Roe flux which satisfies the Rankine-Hugoniot shock jump relations is used 
to calculate the interface fluxes, then the design order of accuracy can be recov- 
ered. Despite its simplicity, this method has serious disadvantages. Detection and 
localization of multidimensional complex shocks, generalization of this approach to 
moving and interacting shocks, and use of ENO type approximations in singular 
regions where the smooth interpolants cannot be constructed make this approach 
unreliable, and, therefore, little used in real numerical applications. 

An alternative strategy is to reduce grid spacing locally near a shock rather 
than refine the grid globally. The idea of this method is based on an assumption 
that the first-order error is generated locally at the shock and then transported 
downstream. This approach, known as a grid adaptation, has been extensively 
used to improve resolution of captured discontinuities for at least twenty years, 
e.g., see [4]. There are two basic strategies of grid adaptation: local grid refinement 
and grid redistribution. In the first approach, grid nodes are added to locally enrich 
the grid to achieve higher accuracy. In the second approach, the number of grid cells 
is fixed and the position of grid points is adjusted to improve the numerical solution 
accuracy. Until now, little attention has been paid to one of the most important 
problems associated with the adaptive grid methods: the essential effect of the grid 
points distribution on error in the numerical solution. It should be emphasized 
that concentration of grid points in regions which most influence the numerical 
solution accuracy may at the same time introduce additional error because of the 
grid nonuniformity [5]. 

Most adaptive grid methods are based on the error equidistribution principle 
developed in [6], [7], [8], which, in turn, is driven by one or another error estimation 
technique. One of the widely used error estimators is gradient or local curvature of 
the numerical solution [9, 10, 11]. An alternative method is to equidistribute the 
local truncation error estimate or the finite element residual which is similar to the 
finite difference truncation error [12, 13]. Another class of error estimators is based 
on evaluating the solution interpolation error [14]. For second-order discretizations, 
this method is reduced to estimation of the local curvature of the numerical solution. 
Richardson extrapolation is also used to estimate error in the numerical solution 
[15]. This procedure compares the solution obtained on the existing grid with one 
computed on a grid that is twice as coarse in each spatial direction. 

Although the error estimators mentioned above are quite different, all of them 
rely on certain smoothness of the differential solution that is not the case for dis- 
continuous flows. In fact, most grid adaptation criteria that can be found in the 
literature become singular at discontinuities. To remove this singularity and to 
make the adaptive grid sufficiently smooth, a grid smoothing procedure must be 
used. As a result, the grid adaptation near discontinuities is driven by the grid 
smoothing procedure rather than the error estimate itself. 
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If grid cell interfaces are not aligned with a shock, the captured discontinuity 
is always smeared over several grid points, which leads to the 0(1) error near the 
shock. Since the numerical solution is first-order accurate away from the shock, the 
true error and its estimates achieve their maximum values at the discontinuity. As a 
result, any grid adaptation procedure based on either the true error or its estimate 
leads to excessive clustering of grid points or local grid refinement near the shock. 
This kind of a grid adaptation is intended to reduce the solution error in the vicinity 
of discontinuities, but it does not necessarily guarantee improvement in accuracy in 
regions where the solution is smooth. As shown in [16], the conventional adaptive 
mesh refinement procedure based on gradient or local curvature of the numerical 
solution can lead to large error in the shock location caused by insufficient accuracy 
in smooth regions of flow ahead of and behind the discontinuity. 

In the present paper, the accuracy of the grid redistribution and local grid refine- 
ment methods are studied. The 2-D test problem used is the supersonic flow around 
a circular cylinder, for which a Chebyshev bow-shock fitting spectral method is em- 
ployed to obtain a very accurate numerical solution [17]. This solution is used as 
the “exact” solution in all subsequent refinement studies. The Euler equations are 
approximated with second- and fourth-order linear shoclc-capturing schemes based 
on the Lax-Friedrichs splitting of the flux vector. The refinement studies show that 
for the second-order scheme, neither grid adaptation strategy improves the numeri- 
cal solution accuracy compared to that calculated on a uniform grid with the same 
number of grid points. For the fourth-order scheme, the dominant first-order error 
component is reduced by the grid adaptation, while the high-order error component 
drastically increases, because of the grid nonuniformity. As a result, the grid adap- 
tation provides improvement in the solution accuracy only asymptotically. A simple 
error analysis of the grid redistribution and local grid refinement methods reveals 
the main reasons why the grid adaptation methods do not improve the accuracy 
of captured discontinuities. Summarizing the numerical and theoretical results, we 
propose a new grid adaptation criterion for captured discontinuities. 

The paper is organized as follows. Section 2 presents the blunt body problem 
and the spectral solution. Section 3 presents the finite difference methods used 
in the studies. Section 4 presents the grid adaptation strategies, including grid 
redistribution and local grid refinement. Section 5 presents the error analysis for 
the effects of grid adaptation on solution accuracy. Section 6 presents the second- 
and fourth-order results, Section 7 presents the new grid adaptation criterion for 
captured shocks, and Section 8 presents the conclusions of the work. 

2. SETTING OF A PROBLEM AND ITS ‘EXACT” SOLUTION 

We consider the steady-state supersonic inviscid gas flow around a circular cylin- 
der. This test problem is desirable because the bow shock is 2-D, rather than simply 
a 1-D shock rotated into a 2-D reference frame. Furthermore, the problem is simple 
enough to obtain an “exact” solution. 

The conservation law form of the 2-D Euler equations written in Cartesian coor- 
dinates (*, y), as follows: 


<9U dF <9G _ . 

~dT + 7h + Ihj ~ °’ 


( 1 ) 
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is used to describe the flowfield. The variables p, u, v, P, and e are the density, x- 
velocity, y-velocity, pressure, and total specific energy, respectively. The governing 
equations are closed with the equation of state for a perfect gas 


P={l-l)P e 


-(n 2 + , 2 ) 


where 7 is the ratio of specific heats, which is assumed to have a constant value of 
1.4. 

To use different grid adaptation techniques, a differentiable one-to-one coordinate 
transformation , 


T — t 

s =Z(x,y) (2) 

v = v(*,y), 


is applied to map a physical domain with curvilinear boundaries onto a unit square. 
Note that the £ and T) coordinates do not depend on time and, therefore, moving 
grids are not considered in the present study. 

The Euler equations in the curvilinear coordinates (2) can be written in conser- 
vation law form as 


<9U <9F <9G _ 

c)r + d(, ^ dr/ 


(3) 


U = ju, F'a= j (6-F I i„G) , G=j (,,,F | ,, y G) 

where the Jacobian of the mapping is given by 


J ~ d(x,y) ~^ Vy ^ yVx ' 


To close the governing equations, boundary conditions should be specified. Be- 
cause of the symmetry along the body centerline, only half of the domain is con- 
sidered. The following boundary conditions are imposed along the symmetry line: 


n=0 — 0) 


du 

dP 

dp 

dr) 

„=o~ d V 

„=o ~~ d V 


(4) 


On the cylinder surface, the no penetration boundary condition, 


uy ( -vx ( \ (=1 = 0, (5) 

is imposed. At the supersonic inflow, all flow quantities are prescribed. The out- 
flow boundary is chosen so that the outflow is fully supersonic and, therefore, no 
boundary conditions are imposed. 
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A spectrally accurate numerical solution to the blunt body problem described 
above is found by using a Chebyshev bow-shock fitting algorithm [17]. The shock 
position always coincides with the inflow boundary along which the Ranlcine- 
Hugoniot relations are used. The physical domain changes continuously as the bow 
shock moves to its steady state position, while the computational domain remains 
unchanged. The Chebyshev collocation method is employed in both the radial and 
circumferential directions to discretize the Euler equations. The equations are then 
marched in time until a steady state solution is reached. Further details on the 
Chebyshev shock-fitting technique can be found in [17]. As has been shown in 
[2], the spectral solution is exact to at least eight significant digits. This solution, 
which is further referred to as the “exact” solution, is spectrally interpolated to a 
sequence of uniformly spaced grids to evaluate error in finite difference solutions. 

3. SECOND- AND FOURTH-ORDER NUMERICAL METHODS 

Second-order fully upwind and fourth-order upwind-biased linear finite difference 
schemes based on the Lax-Friedrichs flux splitting are used to discretize the Euler 
equations. These approximations can be written in a semi-discrete form as 

A = I (D- F + + D+F-) + i (D- G+ + D* G-) , (6) 

where and are linear finite difference operators in £ and 77 , respectively. 

The Lax-Friedrichs fluxes are given by 

F ± = F±|A” ! °^|U 

g± = g±|a;; ,m '|u, ^ 

where |A^ aa: | and |A” lo:r | are the maximum values over the entire domain of the 
contravariant eigenvalues |w| + c and |fi| + c, accordingly. 

The following second- and fourth-order spatial operators, Dt and , are used 
in the present analysis: 


D 2 uj — 2h ( u j - 2 + 3wj) /g\ 

D t w i = + 4 m j + 1 ^ u i+ 2 ). 

°4 u j = + 6«j- 2 - 18«J-1 + 10«J + 3«j+i) , 

D t U J = - 10uj + 18«j+i - 6uy + 2 + u j+3 ), 

where h is a grid spacing either in £ or in ?/. For the second-order scheme ( 8 ), the 
second-order boundary closure is employed at all boundary points. For the fourth- 
order scheme, the third-order stencils used near the boundaries are the optimal 
stencils derived from nearest neighbor information, biased where possible in an 
upwind direction. These second- and fourth-order methods will be referred to as 
“LF-2-2” and “LF-4-3”, respectively. 

All finite difference grids considered in the present study are constructed such 
that a corner grid point coincides with the stagnation point. It has been found 
that this finite difference formulation may lead to numerical instability if all the 
boundary conditions on the symmetry line are imposed weakly through the flux. 
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This instability is caused by the weak formulation of the boundary condition for the 
v component of the velocity vector. Imposing this boundary condition in the strong 
sense eliminates the numerical instability. Numerical calculations have shown that 
the present formulation is as robust as the staggered formulation typically adopted 
by finite difference algorithms near stagnation points. 

On the supersonic inflow, the entire state vector U|t = o is specified. On the 
outflow boundary, the solution is calculated by using high-order fully upwind ap- 
proximations, i.e. , no boundary conditions are imposed, which is consistent with 
the characteristic analysis for the supersonic outflow. On the impermeable wall, 
the no penetration boundary condition is imposed weakly by solving the Riemann 
problem approximately. Note that at steady state, the normal velocity at the wall 
is nonzero, but converges to zero with an order property consistent with overall 
formulation. 

A three stage explicit Runge-Kutta method is used to drive the solution to steady- 
state. To accelerate the convergence, the implicit residual smoothing method pro- 
posed in [18] is employed. 

4. GRID ADAPTATION METHODS 

The main purpose of the present study is to evaluate the influences of grid adapta- 
tion on overall solution accuracy, in the presence of the first-order error component 
resulting from shock capturing. The grid refinement studies presented in [1, 2] 
show that high-order shock-capturing schemes are first-order accurate downstream 
of shocks. Hence, away from the shock, the pointwise error is proportional to the 
grid spacing, and consequently decays like 0(h) as the grid is refined. At the same 
time, if a grid is not aligned with the discontinuity, any shock-capturing numerical 
scheme gives a discrete shock profile with at least one intermediate point in the 
shock. Thus, on sufficiently fine grids, the true pointwise error achieves its global 
maximum value of 0(1) at the discontinuity. To demonstrate this error behavior, 
the centerline pressure error distribution for the M <*, = 3 blunt body problem cal- 
culated using the LF-4-3 scheme is shown in Fig. 1. As one can see in the figure, 
the true pointwise error is singular at the shock located at x = 1.698. This con- 
clusion is based on the property of shock-capturing schemes and does not depend 
on grid spacing in the vicinity of the shock. As a result, any grid adaptation pro- 
cedure based on the error equidistribution principle and the true error will always 
concentrate grid points or refine the grid near discontinuities of the solution. 

It should be emphasized that the error equidistribution principle is not valid 
for problems with strong discontinuities. Actually, this principle is based on the 
minimization of the integral norm of the solution error or its estimate, which, in 
turn, is obtained as the solution of the Euler-Lagrange differential equation which 
cannot be used for discontinuous functions. 

Furthermore, all error estimators found in the literature, such as estimators based 
on the truncation error, the finite element residual, recovery techniques, and the 
extrapolation, require certain smoothness of the solution of the original differential 
problem, which is not a property of discontinuous flows. As a result, these error es- 
timators are singular in regions where the solution is discontinuous. This singularity 
corresponding to the singularity of the true error at the shock leads to excessive 
grid refinement or clustering of grid points around the shock. Asymptotically, any 
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FIG. 1 . Pointwise error for the M 0 0 = 3 blunt body problem obtained with the LF-4-3 
scheme on uniform grids. 


grid generator based on the error equidistribution principle and one of the error 
estimation procedures mentioned above would generate an adaptive grid such that 
the ratio of the local grid spacing at the shock to the neighboring one, generated in 
regions where the solution is smooth, approaches zero. This degeneration in grid 
spacing occurs because these error estimators become singular at the discontinu- 
ity. In practical applications, this singularity is eliminated by excessive smoothing 
of the error estimate function. As a result, the grid adaptation in the vicinity of 
the shock is driven by the smoothing procedure rather than by the error estimate 
itself. Therefore, instead of considering different error estimators and different grid 
smoothing techniques, we generate C 00 adaptive grids which are clustered or lo- 
cally refined near the shock. This kind of a grid adaptation is intended to reduce 
the first-order error component caused by the shock-capturing procedure and to 
increase the overall solution accuracy towards the design accuracy of the numerical 
scheme used. 

Because we know the exact solution, and consequently the exact shock location, 
an adaptive grid can be generated analytically. Quasi-one-dimensional grid adap- 
tation can be employed, assuming that one family of grid lines is aligned with the 
bow shock. It can be done because the numerical solution error in the circumfer- 
ential direction, which is aligned with the shock, is much smaller than that in the 
radial direction. To demonstrate this property, we solve the blunt body problem 
at M O0 = 3 by using the LF-4-3 scheme on tw'o sequences of uniform grids. The 
first one is obtained by doubling the number of grid points in both the radial and 
circumferential directions, as follows: 33 x 33, 65 x 65, 129 x 129, and 257 x 257. 
The second one is generated by refining the grid only in the radial direction, while 
the number of grid points in the circumferential direction remains fixed and equal 
to 33 for all the grids considered. Figure 2 shows that the numerical solution errors 
obtained on these grids are practically identical. It indicates that the discretiza- 
tion error in the direction parallel to the shock is several orders of magnitude less 
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FIG. 2. Grid refinement study at Afoo = 3, showing the pressure error obtained with the 
LF-4-3 scheme on two sequences of grids: 1 ) uniformly refined in £ and 77 , and 2) refined in 
while the number of grid points in rj is fixed and equal to 33. 


than that in the radial direction. Since the numerical solution error is strongly 
dominated by the radial error component, we use 1-D grid adaptation with the 
grid points constrained to move along one family of fixed radial coordinate lines, 
such that the grid is always aligned with the exact bow shock. Note that this grid 
adaptation procedure does not produce skewed cells where the solution accuracy 
may deteriorate because of small values of the Jacobian [4]. 


4.1. Grid Redistribution Method 

Uniform meshes employed for finite difference calculations are constructed by 
using polar coordinates: 


x = — r cos 9 
y = rsin 9. 


( 10 ) 


A mapping from the polar to the computational coordinates is given by 


c — | r -LL 

8 - 1 r o (0)-rl 

V — 9/6 mei x, 


( 11 ) 


where r o, ri, and $ max are the radius of the inflow boundary (which, in the case 
of the exact solution, coincides with the bow shock), the radius of the cylinder, 
and the maximum value of 6, respectively. Note that the inflow boundary has been 
chosen so that the bow shock standoff distance is equal to 3/4(ro($) — ry) for all 
meshes considered. 

As mentioned above, an adaptive grid is constructed by redistributing grid points 
in 1-D fashion along the fixed radial lines. The following coordinate transformation 
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FIG. 3. Ratio of the adaptive grid spacing to the corresponding uniform grid spacing. 


is used to generate the adaptive grid: 


t = 


r i 

f f(r,t)dr 

r 

r l 

f f(r,0)dr 

r o( 6 ) 


V — ^/^max 


( 12 ) 


with the monitor function f(r, 9) defined as 


/M) = - 

d-3 = d.\ + 


/ 


1 + [<^ 3 (^ 2 — 1 ) 

d 2 — 

8 (e)) > 


l + e 5 ( r 


~ r s(e) 

■ s (e)-n 


2 


(13) 


where r s {9) is the bow shock radius, d 1 , d^, and <5 are parameters. Note that the 
same monitor function Eq. (13) is used to generate all adaptive grids, independently 
of the number of grid points. Such a choice of the monitor function provides that 
the metric coefficients are (7°° functions in the entire physical domain, which is 
important for constructing high-order accurate finite difference approximations on 
nonuniform grids. The ratios of the adaptive grid spacing to the uniform one 
corresponding to three sets of the parameters d\ and d% (dj 1 ' = 6.45 and d^ = 2.15, 
d^ = 9.3 and djf' 1 = 3.1, d^ = 12.6 and d[ 3 ^ = 4.2) used in the numerical 
calculations are shown in Fig. 3. The figure shows that the adaptive grid spacing 
achieves its minimum value at the shock, while away from the shock, it behaves 
as a linear function. This linear dependence, which corresponds to the well-known 
exponential stretching, provides that in smooth flow regions, the leading truncation 
error terms resulting from the grid nonuniformity and the approximation of are 
of the same order. 
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4.2. Local Grid Refinement Method 

As shown above, for the blunt body problem in question, the radial error compo- 
nent is much larger than the circumferential one. Therefore, the local grid refine- 
ment. method is also employed in a ID manner. The resulting computational grid 
is obtained as a sequence of structured quadrilateral nested grids that are aligned 
with the shock and arranged in block structures, such that each of them has the 
same family of radial lines, = const. This local grid refinement procedure en- 
sures that there are no “hanging” nodes on the embedded interfaces and that the 
smoothness of the original background grid is maintained. Taking into account that 
the maximum pointwise error occurs at the shock and that the exact shock position 
is known a priori, a family of the nested grids is generated only in the vicinity of 
the exact shock location, while the original uniform grid is used in the rest of the 
domain. The width of the region where the grid is locally refined is an adjustable 
parameter that is chosen to be larger than a stencil of the numerical scheme used. 

This local grid refinement procedure is consistent with the error equidistribution 
principle which says that the grid point distribution is asymptotically optimal if 
some error measure is equally distributed over the field. However, for problems with 
shocks, the error cannot be equidistributed, because the error function is singular at 
the discontinuity. In this case, the stopping criterion for the local grid refinement is 
not trivial. In the present analysis, the adaptive mesh refinement method is used to 
recover the design accuracy of a high-order shock-capturing scheme downstream of 
the shock. This condition can be used as the stopping criterion, which determines 
how many levels of refinement are required to guarantee that the numerical solution 
error obtained on the adaptive grid is consistent with the design order of the scheme 
employed. Using the heuristic model expressing the first-order nature of any high- 
order shock-capturing scheme [2], one can represent the total error in the numerical 
solution as follows: 


e = Cih +C p h p , (14) 

where h is a grid spacing, e is the solution error, p is the design order of the 
numerical algorithm, and Ci and C p are problem dependent constants. If the 
solution smoothness is consistent with the order of the approximation used, then 
C i = 0 and C p is proportional to the (p + l)th derivative of the solution. On the 
other hand, if there are unresolved features in the flow such as shock waves and 
contact discontinuities, the first-order error component is generated by the shock- 
capturing procedure, so that C i ^ 0. Note that the shock error component may 
in general include not only the first-order term, but also higher order terms in its 
expansion. Asymptotically, the shock error component is dominated by the first- 
order term C\h, whereas the design-order error component is dominated by the 
leading truncation error term C p h p . Therefore, only the leading shock error and 
design-order error terms are retained in Eq. (14). 

As follows from Eq. (14), to obtain the pth order accurate solution, the local grid 
spacing near the shock must be of the order of H p , where H is a grid spacing of the 
background mesh. Despite the simplicity of this criterion, the total number of grid 
points, and consequently the complexity of the algorithm, can increase significantly. 
Actually, if a pth-order scheme is used to discretize the governing equations on the 
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background uniform mesh with a grid spacing U, the local grid size in the vicinity 
of the shock should be 


h = 0(H P ) 


(15) 


to obtain a convergence rate of p. 

Assuming that the local grid refinement procedure starts from the background 
mesh and a nested grid at the next level is refined in the radial direction by a factor 
of 2, the total number of levels of refinement needed to satisfy Eq. (15) is 


l = ~{p~ l)log 2 #- (16) 

Upon completion of this local grid refinement procedure, the total number of grid 
points becomes 


N total = l+jj + 2lN„, (17) 

where JV W is a half-width of a region in which the background mesh is locally 
refined. For example, if one uses a fourth-order scheme on a 65 x 65 background 
mesh, the locally refined grid satisfying condition (15) and corresponding to N w = 8 
has 353 grid points in the radial direction. Thus, the total number of grid points 
is increased by a factor of 6. In practical applications, a grid is refined in both £ 
and 77; therefore, the same increase in the number of grid points should be in each 
spatial direction. Note that the time step At is also refined by the same factor 
so that ratios At/ At; and At/ Ai] are the same on all nested grids, ensuring the 
stability with explicit finite difference schemes. As follows from the above example, 
the local grid spacing near the shock would be of the order of 0( 2 -24 ), which makes 
the grid refinement study practically impossible. Therefore, instead of satisfying 
Eq. (15), we perform three series of calculations on locally refined grids that have 
a different number of levels of refinement. During the grid refinement study, the 
background mesh is refined globally. From Eq. (15) it follows that the number 
of levels of refinement is increased by one. For instance, if a locally refined grid 
corresponding to a 33 x 33 background mesh has three levels of refinement, then 
on the next 65 x 65 background grid, one level of refinement is added, so that the 
total number of levels becomes four. 

Remark 1. The grid adaptation procedures described above are based on the 
fact that the exact shock location is known a priori. However, in practice, the 
exact shock position is unknow'n, which may introduce an additional error in the 
numerical solution. As has been shown in [16], traditional grid adaptation methods 
can provide that shocks are well resolved, but their locations are highly inaccurate 
due to the lack of resolution of smooth portions of the solution. 

5. ANALYSIS OF THE ERROR CAUSED BY GRID 
NONUNIFORMITY 

5.1. Error Introduced by Grid Redistribution 

Let us estimate the error in the smooth portion of the numerical solution, which 
is introduced by concentrating grid points near the shock. Because grid points are 
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redistributed in a 1-D manner in the present study, we consider the following 1-D 
scalar equation: 


il>t + fx(il’) = o, (18) 

which can be treated as a 1-D analog of the hyperbolic portion of the 2-D steady 
Euler equations (1). It can be shown by rewriting the Euler equations (1) in non- 
conservative form 


<9U 6>U 

dx +A dy 


= 0 , 


(19) 


where A = (§§■) (§§)• The type of the equations (19) is governed by the 

eigenvalues of the A matrix, which are 


x _ ut/icxAP+tP-c- ( 2 0) 

Am - £2372 , 

where c is the speed of sound. For the blunt body problem under consideration, 
eigenvalues Ai and A 2 are always real, regardless of the local Mach number. There- 
fore, the first-order error generated by the shock-capturing procedure at the shock 
propagates downstream along the streamlines v/u. 

Assuming that the initial profile ip(Q,x) is discontinuous, the solution error of Eq. 
(18) downstream of the captured shock consists of the first-order error component 
caused by the shock capturing and the design-order truncation error component. 
Thus, 


e = CiAx + T P (x), (21) 

where C\ is a constant dependent on the solution, Ax is a local grid spacing, and T p 
is the leading truncation error term of a pth-order shock-capturing scheme used for 
approximation of Eq. (18). Since our aim is to estimate the influence of the spatial 
nonuniformity on the numerical solution accuracy, the temporal error component 
is not considered. To the authors’ knowledge, very little work has been done to 
estimate the coefficient C\ . In the present analysis, it, is assumed that the coefficient 
C\ approaches a constant on sufficiently fine grids. This assumption is corroborated 
by the numerical calculations presented in Section 6. 

To evaluate the leading truncation error term on a nonuniform grid, ap+1 times 
differentiable one-to-one coordinate transformation, 


t = T 
x = *(£), 

between the physical and computational spaces is considered. The nonuniform grid 
in the physical space is generated as images of nodes of a uniform mesh in the 
computational domain Q: 


Xj — x(^j). 


j_ 

J’ 


j = 0, J. 


( 22 ) 
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It is assumed that x(f) is a C°° function such that > 0 £ Q. Transferring 

the ^-derivative in Eq. (18) to the computational space, Eq. (18) is rewritten as 


i/’t + 


jkM 


= 0. 


(23) 


Approximating and x £ by some pi, ii -order finite difference formulas yields 


4 P) (/,) 


j+n 2 

E a nfn 

n=j-n i 


j+m 2 

E 



m=j—m i 


(24) 


where L/, is a finite difference operator; n i,n. 2 , a n and mi, m 2 , fi m depend on 
particular approximations used for evaluating and x^, respectively. It should 
be stressed that the present analysis is performed in regions where the solution is 
smooth, so it is assumed that all derivatives needed for the derivation are continuous 
functions on £ 6 [0, 1]. Expanding the numerator and denominator of Eq. (24) in 
a Taylor series about and omitting the index j on the right hand side, we have 

’Z Clnfn = k I qf/f +1) A^ + 0(A£ P+1 ) 

17L7 1 ( +1) ( 25 ) 

E PmXm = x ( + + 0(A^ P+1 ), 

m=j — m 1 


where = f^pf, + ■ A£ = 1 /./. and Cj, and Cf are constants 

dependent on a n and /3 m , respectively. Substituting Eq. (25) into Eq. (24) and 
taking into account that x ^ >0, 6 [0, 1], one can write 


L h (f s ) 


,4 + qfA^,/f +1) 

* t (i + q?4f*? +1) ) 


+ 0(A£ P+1 ). 


(26) 


If is chosen to be sufficiently small so that A£ p |a:^ ’ +1 */«^| 1. Eq. (26) can 

be linearized as follows: 


L h (f*) = -(fe+CJAZ p ff +1) ) — +0(A£ P+1 ). (27) 

x t v ' \ x t / 

Note that the error introduced by the linearization is of the order of 0(A^ 2p ). 
Neglecting higher order terms in Eq. (27). the leading truncation error term can 

be recast as 


Ap+ i) x (p+ i) 

T p (x) = L h [f x ) -f x = qf Af F — (28) 

X ^ 

With the first-order error component, the solution error downstream of the shock 
becomes 


AP+ 1 ) x ip+V> 

e = CiXt Af + C-p At; p A— — - q*A' p 


( 29 ) 
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From the above equation it follows that the error in the smooth portion of the solu- 
tion consists of three parts. The first one is due to the shock-capturing procedure. 
The second one, which also exists on uniform meshes, arises from the evaluation 
of /f. The third one is caused by the grid nonuniformity. When the grid is clus- 
tered in the vicinity of the discontinuity, the metric coefficient changes rapidly, 
achieving its minimum value at the shock. As has been shown earlier, any grid 
adaptation technique based on the error minimization or the error equidistribution 
principle concentrates grid points near discontinuities of the solution, because the 
true error and its estimates reach their maximum values at the singularity. From 
Eq. (29) it follows that this kind of grid adaptation reduces the first-order error 
component, but at the same time, the second and, especially, the last term of the 
numerical solution error increase drastically in regions where the a:l l>+ 1 becomes 
very large because of the strong grid nonuniformity. In other words, traditional grid 
adaptation methods based on the grid redistribution technique transfer the error 
from the first-order term to the design-order term. Although, asymptotically, the 
solution error is dominated by the first-order shock error, numerical calculations 
show that the coefficient C\ in Eq. (29) may be very small. As will be shown in 
Section 6, when the blunt body problem is solved by using the second-order fully 
upwind scheme based on the Lax-Friedrichs flux splitting, the numerical solution 
error downstream of the shock is dominated by the design-order error component 
for all the grids considered, i.e. , ||Cih|| < ||C p h p ||. As a result, the first-order error 
reduction owing to the grid adaptation is much less than the high-order error in- 
troduced by clustering grid points near the shock, which ultimately decreases the 
overall solution accuracy. 

5.2. Error Introduced by Local Grid Refinement 

In contrast to the grid redistribution method, the local grid refinement approach 
described earlier does not change the metric tensor of the mapping and, there- 
fore, the error component resulting from the high-order derivatives of the metric 
coefficients remains the same as on the nonadaptive background mesh. Despite 
this advantage, the local grid refinement introduces additional error caused by an 
abrupt change in grid spacing at interfaces between coarse and fine grids and by 
interface boundary conditions. 

To illustrate the main problems associated with the local grid refinement method, 
we investigate wave propagation properties when there is a discontinuity in grid 
spacing. A scalar linear equation, 

<F.„ + a<F y = 0, (30) 

is considered as a model problem, where a is a positive constant representing tangent 
of the angle between the freestream flow and the x-axis. In the numerical calcula- 
tions, the parameter a was chosen to be 0.34. The model equation which can be 
treated as a scalar analog of the Euler equations (19) is solved on a unit square 
with the following boundary conditions: <E»(0, j/) = e Iwy and $(x,0) = e ~ Iwax . To 
simplify the analysis, we seek the solution of Eq. (30) in the following form: 

Iu/y 


<F(x, y ) = y?(x)e' 


(31) 
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where I = \/—l. Substituting Eq. (31) in Eq. (30) yields 

ip x + Iauxp = 0. (32) 

The above equation is approximated by using a second-order fully upwind scheme in 
the interior and a two-point first-order upwind approximation at the inflow bound- 
ary. First, we consider a uniform grid with grid spacing h.. In this case, the 
corresponding discrete equations are 

3 i 

-<Pj -2<pj-i + -<pj-2 + Iauhtpj = 0, j = 2,J (33) 


and 


pi — po + Iauihfi = 0. 

Introducing a parameter a — auih , equation (33) becomes 

3 

2 


+ 'Pi ~ '2pj - 1 + T.Pi ~ 2 — 0. 


Substituting the solution of the form 


Pi 


; = C\ J 


into the discrete equation (35) yields 


1 


— T a I A " — 2 A + — — 0 . 


The quadratic equation can easily be solved to give 

li^i-Hf + a/) 


Ai o — 


| + al 


Thus, the general solution of Eq. (35) is given by 

Pi = ciAj + c 2 A J 2- 


(34) 


(35) 


(36) 


(37) 


(38) 


(39) 


The unknown coefficients c 1 and c 2 are found by using the initial condition p\ x= o = 
po and Eq. (34), as follows: 


Cl + C 2 = p 0 

(1 + Q'/)(c 1 Ai + c 2 A 2 ) — ci — c 2 = 0. 


( 40 ) 


The exact discrete solution Eqs. (31, 38-40) has been derived under the assump- 
tion that the computational grid is uniform. However, the same approach can be 
applied to construct the discrete solution on a locally refined grid. Without loss of 
generality, only a two-level locally refined grid with a factor of two refinement, i.e., 
h\/h 2 = 2, is considered. The exact discrete solution on the locally refined grid can 
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be constructed as follows. Assuming that the coarse grid with grid spacing h i = h 
covers the left half of the domain, and that the discontinuity in mesh size is located 
at the midpoint x = 1/2, the solution at this point is 

9%=ch\ J A 2 + ci 2 Af 2 /2 , (41) 

where cn = c\, C 12 = C 2 , and An and A 12 are roots of the characteristic polynomial 
Eq. (37) corresponding to the background mesh h,\ — h. The above numerical 
solution can be treated as an incident wave traveling from left to right, which 
is then transmitted by the interface between the coarse and fine grids. For the 
hyperbolic equation (30), the solution of the Riemann problem at the interface is 
trivial fright = y’left- Since the same fully upwind scheme Eq. (33) is used in the 
right half of the domain, the discrete solution on the fine nested grid can be written 
in the same form as Eq. (38-40). However, instead of the grid spacing hi = h, 
one should use the grid spacing h .2 = h/ 2, and the Direchlet boundary condition 
at x = 0 should be replaced with the solution of the Riemann problem at the 
midpoint, i.e. , u| a , = x /2 = y?j/ 2 , where fj/2 is given by Eq. (41). Thus, the exact 
discrete solution on the locally refined mesh is 

'P'j ^ = C21^21 + c 22^22- (42) 

where A 21 and A 22 . are defined by Eq. (38) with a = awho, and the coefficients C 21 
and C 22 are the solution of the following linear system of equations: 

c 2 i + c 22 = fjj 2 ( 43 ) 

(1 + Cv/) ( O21 Ao 1 + C22A22) — C 21 — c 22 = 0 . 


The exact discrete solution Eqs. (42-43) calculated on the locally refined grid is 
compared with the uniform grid solution at x = 1. To determine the influence 
of the grid discontinuity on the numerical solution accuracy, we compare the error 
obtained at x = 1 on both the locally refined and corresponding uniform grids. The 
error is calculated as the absolute value of the difference between the discrete and 
analytical solutions of the original differential equation (30), which is 

$(x,y) =e~ Iawx e Iwy . (44) 

The ratio of the error obtained on the locally refined mesh to the uniform 
grid solution error corresponding to four grid spacings of the background mesh 
(h = 10~ 2 , 10~ 3 , 10~ 4 , and 10 -5 ) is plotted versus u i 3 h 2 in Fig. 4. Note that the 
parameter ui 3 h 2 has been chosen so that it is proportional to the leading truncation 
error term of the second-order scheme used. Therefore, for small values of lo 3 h 2 
corresponding to the fully resolved solution, the error ratios calculated on grids 
with different background grid spacings are practically identical. The curves cor- 
responding to different mesh densities begin to diverge at a critical value of u 3 h 2 
that corresponds to the point where the error obtained on the locally refined mesh 
becomes equal to the error obtained on the background mesh. Since this point 
is asymptotically independent of the grid spacing, we can introduce the critical 
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FIG. 4. Ratio of the error obtained on the two-level locally refined grid to the error obtained 
on the corresponding background grid. 


frequency 


(Hf 


(45) 


If 0 < S < u> < w„, where S is a small positive constant, then the locally refined grid 
provides better accuracy than the corresponding background mesh. Note that the 
constant in the formula (45) for the critical frequency u;* depends on the numerical 
scheme used, the angle between the freestream flow and the *-axis, and the size 
of the domain. Despite this dependence, for any stable second-order scheme, the 
qualitative behavior w* ~ h ~ 2 / 3 always remains the same. 

If the frequency of the exact solution exceeds the critical value , it cannot be 
resolved on a uniform grid with the grid spacing h. As a result, the errors ratio 
exhibits oscillatory behavior. Note, however, that for all the grids considered, there 
exists the local maximum at ui 3 h 2 ps 15, corresponding to the frequency that is most 
strongly amplified by the abrupt change in grid spacing, as is evident in Fig. 4. 
The main reason for such a behavior is the hyperbolic nature of Eq. (30) . Actually, 
if error has been introduced on the coarse grid, it propagates downstream along the 
characteristic. In the case of discontinuous solutions, such as shock waves, the cap- 
tured discontinuity involves all frequencies and, therefore, cannot be fully resolved 
on any grid. The error component corresponding to these unresolved frequencies 
propagates downstream along the characteristics and amplifies at interfaces. The 
error amplification occurs because of the discontinuity in grid spacing and the ap- 
proximate interface boundary conditions. It should be emphasized that the error 
amplification becomes stronger, if either the number of interfaces or the grid refine- 
ment factor increases, or the grid is globally refined. 

As one can see in Fig. 4, there is a boundary layer at ui 3 h 2 — >■ 0. This deteriora- 
tion in accuracy for the lowest frequencies is caused by reducing by one the order of 
approximation at the first, point of the coarse and fine grids. To demonstrate this 
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property, we estimate the ratio of the solution error obtained on the locally refined 
mesh to one obtained on the corresponding uniform grid, when the exact solution 
frequency is sufficiently small, i.e., 

w <C /t~ 2/3 . (46) 


If inequality (46) holds, then the truncation error analysis can be applied. Taking 
into account the fact that the error of integration of Eq. (32) is proportional to the 
length of the integration interval, the ratio of the truncation errors can be estimated 
as follows: 


£2 

ft 


[Cw„h 2 + C 2 <p xxx h 2 (| - h)] + [c^ xx (|) 2 + C 2 <p xxx (|) 2 (I 
Ci<p xl! h 2 + C 2 ip xxx h 2 (1 -h) 



(47) 


where C\ and C 2 are constants dependent on the first- and second-order approxi- 
mations used, e.g., in our case C\ = —1/2 and C 2 = —1/3. In Eq. (47), the first 
term in the square brackets is the integration error on the coarse grid with grid 
spacing h, and the second term is the integration error obtained on the fine grid 
with grid spacing h/2. Note that the C\<p xx h 2 and Ci<p xx (h/2) 2 terms are due to 
the first-order approximation used at the first grid points of the coarse and fine 
grids, respectively. With the exact solution Eq. (44), Eq. (47) is reduced to 


€2 = fCl + IC 2 u> (1 - h) + IC 2 % (I - I) 
€i Cl + IC 2 oj (1 — h) 

Letting h — > 0 yields 


e 2 _ 5 

77 - 4 (ci cy) 


\Jct + \c vh \cic w. 


(48) 


(49) 


From the above equation it follows that for the lowest frequency w — t 0, the trun- 
cation error calculated on the globally uniform grid is 5/4 times less than that ob- 
tained on the locally refined grid having one and a half times as many grid cells as 
the corresponding uniform mesh. As the wave frequency increases, the errors ratio 
tends to its asymptotic value of 5/8. Note that if the governing equation is approx- 
imated with a second-order scheme in the entire computational domain, then the 
first-order terms in Eq. (47) vanish, providing that e 2 /ei = 5/8, Vu : u> h~ 2 / 3 . 

Remark 2. In solving the 2-D Euler equations, two additional sources of errors, 
which are not present in the test example considered above, can arise. The first one 
is due to an approximate solution of the Riemann problem at interfaces between 
the coarse and fine grids. The second source of error is caused by the reflection 
of waves going through interfaces. Note that in the above example, there is no 
numerical reflection, because the fully upwind formulation, which is consistent with 
the characteristic of Eq. (30), is employed. 


6. RESULTS AND DISCUSSION 

The inviscid M m = 3 flow' around a circular cylinder is used to test the accuracy of 
the grid redistribution and local grid refinement methods. A grid refinement study 
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is performed to investigate the influence of the grid adaptation on the numerical 
solution accuracy downstream of the captured bow shock. The following sequence 
of grids is used in the grid refinement study: 33 x 33, 65 X 65, 129 x 129, and 
257 x 257. Error in the smooth portion of the solution is measured in the L 2 sense, 
as follows: 


II*-*' 



(50) 


where <pkj and <j >\ :*■ are the numerical and exact values of (j> at point ( k,j ), and K 
and J are the number of grid points in the domain of interest. The norm of wall 
quantities is formed in a manner similar to Eq. (50), but only the wall points are 
used in the formula. 


6.1. Grid Redistribution Method 

An adaptive grid is generated by equidistributing the monitor function Eq. (13), 
as described in Section 4.1. Three sets of the parameters di and do (6.45 and 2.15, 
7.95 and 2.65, 12.6 and 4.2) in Eq. (13) generate three families of the adaptive grids. 
These three sets of parameters d\ and do provide that the ratio of the adaptive 
grid spacing to the uniform grid spacing at the shock is equal to 1.5, 2.5, and 4, 
respectively. A sample 65 x 65 adaptive grid generated by this grid redistribution 
procedure is shown in Fig. 5. The region around the shock is well resolved by 
reducing the local grid spacing by a factor of 2.5 compared with the uniform mesh 
with the same number of grid points. Figure 6 shows the pressure contours obtained 
with the LF-2-2 scheme on the 65 x 65 uniform (left) and adaptive (right) grids. 
As one might expect, the grid adaptation improves the shock resolution. To give 
greater insight into how the grid clustering influences on the shock smearing, the 
centerline pressure distributions computed on these adaptive and uniform grids 



X 

FIG. 5. Adaptive 65 X 65 grid (Aa7 un /Aa: m i u = 2.5) generated by the grid redistribution 
method. 
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X X 

FIG. 6. Pressure contours obtained with the LF-2-2 scheme on 65 X 65 uniform (left) and 
adaptive (right) grids. 



FIG. 7. Centerline pressure distributions calculated using the LF-2-2 scheme on the 65 X 65 
uniform and adaptive grids. 


are depicted in Fig. 7. For all the grids considered, the bow shock is captured 
within 6 cells, which is typical for the Lax- Friedrichs flux splitting employed. At 
the same time, the shock thickness, which is proportional to the local grid spacing 
in the physical space, becomes thinner as the clustering of grid points increases. 
Locally near the shock, the numerical solution profile, which remains practically 
unchanged, is scaled according to the local grid spacing. Note that the numerical 
solutions obtained with the LF-4-3 finite difference scheme demonstrate similar 
behavior. 

Grid refinement studies using three families of adaptive grids and one family of 
uniform grids that have the same number of grid points are presented in Figs. 8-13. 
The Ln norm of the wall pressure error is shown in Fig. 8. As one can see in the 
figure, the convergence rate obtained for the smooth problem is consistent with the 
design order of the scheme used. The smooth problem is formulated in the region 
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FIG. 8. Grid refinement study at Moo = 3, showing the wall pressure error obtained with 
the LF-2-2 scheme on uniform and adaptive grids. 



FIG. 9. Crid refinement study at M 0 0 = 3, showing the pressure error in the half of the 
domain closest to the body obtained with the LF-2-2 scheme on uniform and adaptive grids. 


bounded by the exact shock wave and the cylinder, so that there is no discontinuity 
in the domain. The inflow boundary condition at the bow shock is implemented by 
solving the Riemann problem between the numerical state and the exact post-shock 
conditions from the spectral solution. 

Although the shock resolution is improved by clustering grid points around the 
shock, the error in the pressure on the body surface increases by a factor of 5-10 
compared with that calculated on the uniform grid with the same number of grid 
points. It should be noted that the higher the concentration of grid points in the 
vicinity of the shock, the larger the error that is introduced into the numerical 
solution by the grid nonuniformity. The grid adaptation results in the metric coef- 
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FIG. 10. Grid refinement study at Moo = 3, showing the wall pressure error obtained with 
the LF-4-3 scheme on uniform and adaptive grids. 



FIG. 11. Grid refinement study at Moo = 3, showing the pressure error in the half of the 
domain closest to the body obtained with the LF-4-3 scheme on uniform and adaptive grids. 


ficients drastically increasing near the shock, which, in turn, considerably decreases 
the accuracy in regions where the solution is smooth. These numerical calculations 
corroborate the analysis presented in Section 5.1. Surprisingly, the shock-capturing 
solutions obtained on uniform grids not only are more accurate than the corre- 
sponding solution of the smooth problem, but also exhibit higher convergence rate 
on fine meshes. This anomalous behavior present in the shock aligned cases did not 
generalize to non-shock aligned meshes. Nevertheless, this superconvergence is not 
yet well understood and requires additional investigation. 

The L 2 norm of the pressure error calculated in the half of the domain closest to 
the body is presented in Fig. 9. Such a choice of the domain enables us to exclude 
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any point close to the captured shock, where the point. wise error is of the order 
of 0(1). The error behavior obtained in the field is very similar to that exhibited 
by the wall pressure error. As in the foregoing comparison, the grid adaptation 
significantly increases the error in the smooth portion of the solution. Note that for 
the LF-2-2 scheme, the first-order error component is comparable with the design- 
order error, giving the appearance of the second-order convergence rate on all the 
grids considered. 

Error convergence plots obtained with the fourth-order upwind-biased LF-4-3 
scheme are depicted in Figs. 10 and 11 which are analogous to Figs. 8 and 9, 
accordingly. In contrast to the second-order scheme, the fourth-order method using 
uniform grids quickly approaches the asymptotic limit and exhibits just the first- 
order convergence on fine grids. In spite of the fact that the solution error is 
dominated by the first-order error component, the grid adaptation improves the 
numerical solution accuracy only asymptotically. On one hand, the concentration 
of grid points near the shock reduces the first-order error component, so that the 
convergence rate increases towards the design order limit, as shown in Figs. 10 
and 11. On the other hand, the high-order error component significantly increases 
because of the grid nonuniformity. As a result, the solution error is dominated 
by the design-order error component. Despite some improvement in accuracy on 
the coarsest adaptive grids, the Lo norm of the solution error obtained on finer 
adaptive grids is larger than that calculated on the corresponding uniform grids. 
As shown in Section 5.1, the main reason for such a behavior is that the design- 
order error component, resulted from the approximation of the metric coefficients, 
drastically increases, owing to the clustering of grid points near the shock. Although 
the adaptive grid solution exhibits a higher convergence rate, the errors obtained 
on the corresponding adaptive and uniform grids become comparable only on the 
finest 257 x 257 mesh. Hence, the grid redistribution method based on the LF- 
4-3 scheme provides improvement in accuracy only on extremely fine grids which 
are very seldom found in practical applications, because of the computational cost 
involved. Note that the adaptive grid solution error suffers from the first-order 
degeneration but with much lower first-order error constant, as follows from Eq. 
(29). 

In the grid convergence studies presented above, the solution error is measured 
in regions where the solution is smooth. However, it appears intuitively that the 
grid adaptation should reduce the solution error in the vicinity of the shock. One 
can argue that calculation of the lift or the drag for discontinuous flows results in 
an error in these integral quantities of 0(AP s Ax s ) , where AP. and A* s are the 
pressure jump across the shock and the grid spacing at the shock, respectively. 
This assumption is the basis for using grid adaptation methods to reduce the first- 
order shock error in the lift and drag which are two of the most important integral 
quantities obtained from the solution of the Euler equations. 

To check whether the grid adaptation improves the accuracy of the integral quan- 
tities, we compare errors in the pressure integral across the bow shock along line 
y = 0 calculated on the same adaptive and uniform grids. Two integration al- 
gorithms have been used to compute the integral. The first one is a standard 
two-point second-order trapezoidal rule quadrature formula. The second one uses 
a piecewise cubic spline to fit. the discrete pressure. Then, the spline is integrated 
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FIG. 12. Grid refinement study for the pressure integral error obtained with the LF-2-2 
scheme on uniform and adaptive grids. 



FIG. 13. Grid refinement study for the pressure integral error obtained with the LF-4-3 
scheme on uniform and adaptive grids. 


analytically, providing a fourth-order accurate integration formula for sufficiently 
smooth functions. Although the fourth-order integration procedure imposes more 
severe constraints on smoothness of the integrated function, pressure integral er- 
rors obtained with second- and fourth-order quadrature formulas are practically 
identical. Therefore, only the results calculated using the fourth-order integration 
procedure are presented. 

Figures 12 and 13 show the grid refinement study for the pressure integral ob- 
tained with the LF-2-2 and LF-4-3 shock-capturing schemes, respectively. In con- 
trast to conventional wisdom, the pressure integral calculated with the LF-2-2 
scheme on clustered grids is less accurate than that computed on uniform grids 
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with the same number of grid points. Figure 12 shows that the accuracy of the 
pressure integral deteriorates as the concentration of grid points in the vicinity of 
the shock increases. For the LF-4-3 scheme, slight improvement of the pressure 
integral accuracy can be observed on both the coarsest and finest adaptive grids, 
while similar deterioration in accuracy occurs on 65 x 65 and 129 x 129 meshes. 

It should be emphasized that the error in the pressure integral across the shock 
behaves very similarly to the L 2 norm of the wall pressure error. This is no surprise 
because these quantities are closely connected. To show this relation, we integrate 
the 2-D steady Euler equations over the entire domain to give 

J (Fj, + G y ) dV = /V ndfi = 0, (51) 

v n 


where T is the flux tensor, fi is the boundary of the physical domain V, n is an 
outward pointing unit vector normal to ft. As follows from Fig. 5, ft = AB + BC + 
CD + DA. Thus, 


j T ndft = 

AB 


J T ■ nd.il 

BC 


j T ndft 

CD 



DA 


T ■ ndft. 


(52) 


Because the finite difference schemes used are fully conservative, Eq. (52) holds not 
only for the exact analytical solution, but also for the discrete solution. As follows 
from Eq. (52), the accuracy of the solution integral along the centerline AB is 
determined by the accuracy of the integrals calculated along the other boundaries 
including the body surface BC. This indicates that the integral Ln norm of the 
wall pressure error is closely related to the accuracy of the pressure integral along 
the centerline. 


6.2. Local Grid Refinement Method 

As follows from the analysis presented in Section 5.2, the most troublesome parts 
of the local grid refinement method are the grid discontinuity and the interface 
boundary conditions. Therefore, special attention has been paid to implementation 
of the interface boundary conditions and their accuracy. In the present study, 
the Roe’s approximate Riemann solver is used to build a flux at grid interfaces. 
The left and right states at a grid interface are obtained as solutions of the Euler 
equations calculated separately in each subdomain. The numerical flux built this 
way is then used to construct the second- and fourth-order approximations, Eqs. (8- 
9). This implementation of the interface boundary conditions maintains stability, 
conservation, and accuracy in multiple dimensions for both second- and fourth- 
order methods. Because only matching nested grids are considered, there are no 
hanging nodes in the domain. 

To verify the accuracy of the interface boundary' conditions described above, two 
calculations of the smooth M co = 2.5 blunt-body problem are performed. The first 
grid refinement study is done on a sequence of uniformly spaced grids. Each grid 
is divided on eight subdomains containing the same number of grid points in both 
the radial and circumferential directions, as shown in Fig. 14. Pressure contours 
computed with the LF-4-3 scheme on this multiblock grid are also shown in Fig. 14. 
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FIG. 14. 8-block 25 X 33 uniform grid (left) and pressure contours (right) of the smooth 
Moo — 2.5 flow calculated with the LF-4-3 scheme. 



FIG. 15. Two-level 41 X 33 locally refined grid (left) and pressure contours (right) of the 
smooth Mqo = 2.5 flow calculated with the LF-4-3 scheme. 


Although the present formulation of the interface boundary conditions results in 
the numerical solution being discontinuous at interfaces, the absolute value of the 
solution jump at the interface is of the order of the truncation error of the scheme 
used. For example, the characteristic interface pressure jump calculated with the 
LF-4-3 scheme on a 25 x 33 uniform 8-block grid is O(10 -3 ). 

The second series of calculations is performed on a sequence of locally refined 
grids similar to those which are used in the shock-capturing formulation. As in 
the previous test case, only the flow between the bow shock and the cylinder is 
considered, i.e., there are no discontinuities in this subproblem. A three-level locally 
refined grid corresponding to a 25 x 33 background mesh and the pressure contours 
calculated on this grid are shown in Fig. 15. Qualitatively, the solutions obtained on 
uniform and nested multiblock grids are very similar. However, the corresponding 
L 2 pressure error norms calculated in the half of the domain closest to the body 
are quite different, as seen in Fig. 16. As shown in Section 5.2, the presence of 
interfaces in the domain introduces additional error in the numerical solution. On 
the coarsest 8-block uniform grid, the solution is less accurate than that obtained on 
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FIG. 16. Crid refinement study at Mqc = 2.5, showing the pressure error in the half of 
the domain closest to the body obtained with the LF-4-3 scheme on the 8-block uniform and 2- 
to 5-level locally refined grids. 


the corresponding single-block uniform grid. This reduction in accuracy is due to 
the fact that all grid points in each subdomain are treated as the boundary points. 
Therefore, the third-order boundary closure approximation is used on the entire 
coarsest mesh. On finer meshes, the multiblock uniform grid formulation exhibits 
the design-order convergence rate and provides practically the same accuracy as 
on the single-block uniform grid. In contrast to multiblock uniform grid results, 
the locally refined grid formulation leads to both one order of magnitude reduction 
in accuracy and deterioration of the convergence rate to 3 on fine meshes. These 
numerical results corroborate the theoretical analysis which shows that the error 
components corresponding to the lowest fully resolved and the high unresolved 
frequencies are amplified by the grid discontinuity. Such a reduction in accuracy is 
also caused by the approximate solution of the R.iemann problem and by spurious 
reflection of waves traveling through the interfaces. 

To test the accuracy of the local grid refinement method, the LF-2-2 and LF-4-3 
schemes are used to capture the bow shock around a Mach 3 circular cylinder. A 
typical three- level 137 x 65 locally refined grid used in both second- and fourth- 
order formulations is shown in Fig. 17. This adaptive grid corresponding to 65 x 65 
background uniform mesh is locally refined near the exact shock location, such 
that the local grid spacing of the finest nested subgrid is 8 times smaller than 
that of the background mesh. Figure 18 shows the pressure contours calculated on 
these uniform and locally refined grids. As expected, the shock w'ave front becomes 
much sharper on the adaptive grid. To demonstrate the improvement in the shock 
resolution owing to the local grid refinement, the centerline pressure distributions 
computed on the three-level and five-level locally refined grids are compared with 
the uniform grid results in Fig. 19. The comparison shows that the bow shock is 
smeared over six grid cells for all the grids considered, while the shock thickness in 
the physical space is reduced proportionally to the local grid spacing. 
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FIG. 17. Three-level 137 X 65 locally refined grid. 




FIG. 18. Pressure contours obtained with the LF-4-3 scheme on the 65 X 65 uniform (left) 
and 3-level 137 X 65 locally refined (right) grids. 


Despite improved shock sharpness, there is no reduction in the L 2 error norm 
measured in smooth portion of the numerical solution calculated with the second- 
order method, as is evident in Figs. 20 and 21 which are analogous to Figs. 8 and 9. 
Figure 20 shows that the L 2 norm of the wall pressure error obtained with the LF- 
2-2 scheme on locally refined grids increases by a factor of 10 1 — 10 2 compared with 
the uniform grid results. The pressure error measured in the half of the domain 
closest to the body exhibits similar behavior, as is evident in Fig. 21. Note that 
for the second-order scheme, the error convergence rate of the local grid refinement 
method is less than the design order. 

In contrast to the second-order method, the numerical solution error obtained 
with the LF-4-3 scheme is dominated by the first-order error component. Figures 
22 and 23 show that the local grid refinement near the shock reduces the first-order 
error component generated by the LF-4-3 shock-capturing procedure, so that the 
design-order error component becomes dominant on the adaptive grids. Although 
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FIG. 19. Pressure distributions along line y = 0 calculated on the 65 X 65 uniform and 3-, 
and 5-level locally refined grids. 



FIG. 20. Crid refinement study at A/qo = 3, showing the wall pressure error obtained with 
the LF-2-2 scheme on uniform and locally refined grids. 


the error convergence rate is almost recovered to the design order, the L 2 norm of 
the wall pressure error obtained on locally refined grids becomes comparable with 
the uniform grid results only on the finest mesh. For coarser grids, the superiority 
of uniformly spaced grids is evident. The same conclusion can be drawn for the 
pressure error measured in the field away from the shock, except that the solution 
errors obtained on the locally refined and uniform grids become comparable on 
the coarser 129 X 129 grid, as one can see in Fig. 23. The numerical results 
obtained with both the second- and fourth-order schemes show that as the number 
of interfaces between coarse and fine meshes increases, the deterioration in accuracy 
also increases, which qualitatively corroborates the analysis presented in Section 5.2. 
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-■ Uniform grid, shocked problem 

-• Uniform grid, smooth problem 



FIG. 21. Grid refinement study at M 0 0 = 3, showing the pressure error in the half of the 
domain closest to the body obtained with the LF-2-2 scheme on uniform and locally refined grids. 



FIG. 22. Grid refinement study at A/qo = 3, showing the wall pressure error obtained with 
the LF-4-3 scheme on uniform and locally refined grids. 


Error convergence plots of the pressure integral along the centerline calculated 
using the LF-2-2 and LF-4-3 schemes on the same sequences of locally refined 
and uniform meshes are depicted in Figs. 24 and 25, respectively. Note that the 
error convergence on uniform grids is monotonic, whereas the error convergence on 
adaptive grids exhibits nonmonotonic behavior. Figure 24 shows that the local grid 
refinement method based on the LF-2-2 scheme does not reduce the pressure integral 
error compared with the uniform grid results. Certain improvement in accuracy of 
the pressure integral computed with the LF-4-3 scheme can be observed on fine 
locally refined grids, as shown in Fig. 25. However, this gain in accuracy decreases 
as the background mesh is globally refined, and it practically vanishes on the finest 
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FIG. 23. Grid refinement study at M 0 0 = 3, showing the pressure error in the half of the 
domain closest to the body obtained with the LF-4-3 scheme on uniform and locally refined grids. 


Uniform grid, shocked problem 



FIG. 24. Grid refinement study for the pressure integral error obtained with the LF-2-2 
scheme on uniform and locally refined grids. 


adaptive mesh. Although only the error in the pressure and its integral have been 
presented in this study, similar convergence behavior is observed for the other flow 
quantities. 

The results presented above have been obtained under the assumption that one 
family of grid lines is parallel to the bow shock. However, in practical applications, 
generation of shock-aligned grids is a complicated problem, because the exact shock 
location is unknown. To gain greater insight into the error behavior when the grid 
is not aligned with the shock, we perform a grid refinement study on uniform grids 
generated by the mapping Eq.(ll) with ro((?) = const. In this case, circumferential 
grid lines are circular arcs which are not aligned with the bow shock. Figure 26 
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FIG. 25. Grid refinement study for the pressure integral error obtained with the LF-4-3 
scheme on uniform and locally refined grids. 



shows a uniform 65 x 65 grid and isobars of the blunt body flow calculated with the 
LF-2-2 scheme on this grid. Note that the maximum angle between the bow shock 
and circumferential grid lines, which occurs at the point where the shock crosses 
the outflow boundary, is about 45°. The L% norm of the wall pressure error and the 
pressure error norm measured in the half of the domain closest to the body surface 
(obtained on the nonaligned uniform grids) are compared with the shock-aligned 
uniform and most accurate adaptive grid results in Figs. 27 and 28, respectively. 
As follows from the grid refinement study, when the grid is aligned with the bow 
shock, the coefficient in front of the first-order error component is negligibly small, 
giving the appearance of a second-order convergence rate. However, this coefficient 
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-u Aligned uniform grid, shocked problem 

-# Aligned uniform grid, smooth problem 



FIG. 27. Grid refinement study at Moo = 3, showing the wall pressure error obtained with 
the LF-2-2 scheme on nonaligned uniform, shock-aligned uniform and adaptive grids. 


-■ Aligned uniform grid, shocked problem 

-# Aligned uniform grid, smooth problem 



FIG. 28. Grid refinement study at Moo = 3, showing the pressure error in the half of the 
domain closest to the body obtained with the LF-2-2 scheme on nonaligned uniform, shock-aligned 
uniform and adaptive grids. 


becomes much larger on nonaligned grids that makes the first-order shock error 
component dominant on much coarser grids. Despite the fact that the L 2 norm of 
the pressure error obtained on the finest nonaligned uniform mesh is one order of 
magnitude less accurate than that on the corresponding shock-aligned uniform grid, 
the nonaligned uniform grid solution is still more accurate than the best adaptive 
grid results. As one can see in Figs. 27 and 28. the error convergence obtained 
with the LF-2-2 scheme is very similar to that calculated with the LF-4-3 scheme 
on shock-aligned grids. Figure 28 shows that the grid adaptation can improve the 
numerical solution accuracy only asymptotically. 
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7. GRID ADAPTATION CRITERION FOR CAPTURED SHOCKS 

Based on the results presented in Sections 6.1 and 6.2, we propose the following 
grid adaptation strategy for captured shocks. The most general error convergence 
behavior shown schematically in Fig. 29 is characterized by the presence of three 
different regions in the error convergence plot. The first, region corresponds to 
coarse meshes such that the shock error component is of the order of 0(1). In 
this region, the grid adaptation near the shock is desirable because it reduces the 
most troublesome part of the error generated by the shock-capturing procedure. It 
should be stressed that the grid adaptation not only reduces the shock error, but 
also localizes the captured shock, which minimizes interaction of Gibbs oscillations 
with the boundary conditions and other features of the flow. This is one of the main 
reasons why high resolution shock-capturing schemes, such as ENO and weighted 
ENO schemes, provide better accuracy on coarse grids compared to high-order 
linear algorithms. 

In the second region, which corresponds to sufficiently fine grids, the coefficient 
C\ in Eq. (14) approaches its asymptotic value which, as has been shown for 
the LF-2-2 scheme, may be very small. For these grids, the design-order error 
component dominates the first-order error component, i.e. ||C p /? l: ’|| > |jCi/i||. As 
a result, any clustering of grid points or local mesh refinement in the vicinity of 
the shock reduces only that part of the error which is not dominant, ultimately 
reducing the numerical solution accuracy in smooth portions of the solution. The 
main reason for such a behavior is the fact that both grid adaptation strategies 
are not without penalties in the solution accuracy. As follows from the truncation 
error analysis, for the grid redistribution method, the main source of error is due 
to the high-order derivatives of the metric coefficients which drastically increase on 
nonuniform meshes. The numerical analysis performed for the local grid refinement 
technique has shown that, in this case, the error accumulation occurs because of 
discontinuities in grid spacing and the coupling between coarse and fine meshes, 
herein calculated as an approximate solution of the Riemann problem. Because 



FIG. 29. Error convergence diagram for a high-order shock-capturing method. 
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of these penalties in the solution accuracy, the grid adaptation increases the total 
error in the numerical solution if it is used in the second region. 

The third region corresponds to the first-order asymptotic limit. In this region, 
the grid adaptation near the bow-shock would be the most efficient use of the com- 
puter resources. However, as has been shown numerically, there is a transition zone 
where the error convergence rate is 0(h ), but the grid adaptation does not provide 
improvement in the numerical solution accuracy compared with the corresponding 
uniform grid results. This is because both the grid redistribution and local grid 
refinement methods introduce additional error in the numerical solution, caused by 
the grid nonuniformity. Note that the transition zone becomes larger if the com- 
putational cost is used instead of the number of grid cells in the error convergence 
plots. 

From the above analysis it follows that the grid adaptation is desirable if the 
following inequality holds: 


JICiAH » \\C p h p \\, (53) 

where the norm is measured in regions where the solution is smooth. The above 
inequality can be used as a grid adaptation criterion for captured discontinuities. 
It should be noted that the constants C\ and C p depend on both the problem 
and the numerical scheme used. As has been mentioned earlier, the coefficient 
C p can be treated as the leading truncation error term. To our knowledge, there 
are no theoretical results for evaluation of the coefficient C 1 . Therefore, to use 
the criterion (53) in practical applications, global grid refinement or coarsening is 
required. The global coarsening can be used if the grid is fine enough to correspond 
to the second and third regions in the error convergence diagram. However, if the 
grid is very coarse that corresponds to the first region, the coarsening cannot be 
used, and the grid refinement is tile only way to evaluate the error convergence 
rate. This approach becomes quite expensive in three dimensions. From this point 
of view, it is very important to be able to predict a priori the error introduced by 
the shock-capturing procedure. This quantification together with the criterion (53) 
will provide guidance for grid adaptation for captured discontinuities. 

It is instructive to speculate on the generality of this work, particularly on two 
practical questions: 1) are the conclusions and suggestions (see Fig. 29) valid for any 
general discontinuous flow, and 2) are the present results (second- and fourth-order 
Lax- Friedrichs schemes) general for any high-order numerical algorithm? Although 
a precise answer does not presently exist to either question, we make the following 
conjecture: Similar qualitative results and an equivalent grid adaptation criterion 
(53) will be obtained for any flows containing multidimensional shocks, indepen- 
dent of the high-order numerical method (central difference with scalar dissipation, 
Roe, TVD, and ENO). This assertion is based on the observation that for any 
captured discontinuities, the numerical solution error obtained with any high-order 
shock-capturing method consists of two parts: the first-order shock error C\h and 
the design-order error component C p h F . Although the coefficients C\ and C p are 
problem and scheme dependent, asymptotically, these coefficients do not depend 
on the grid spacing h. Consequently, for any high-order method, such that p > 1, 
and any finite values Ci and C v , the solution error is asymptotically dominated by 
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the first-order error component, corresponding to the third region in Fig. 29. If 
C i <StC C p , the solution error exhibits high-order convergence on sufficiently coarse 
grids, and corresponds to the second region in the error convergence diagram. If 
the grid is very coarse, such that C\h ~ 0(1), the solution error corresponds to 
the first region in Fig. 29. Note, however, that the first and second regions may 
be very small or even nonexistent if (\ C p , whereas the third region always ex- 
ists, regardless of the discontinuity strength and the high-order numerical algorithm 
used. 

Although the results presented above have been obtained for steady state flows, 
the same conclusions can be drawn for time-dependent problems as well. An essen- 
tial effect of the grid nonuniformity on the numerical solution accuracy remains one 
of the most important sources of error for unsteady problems with shocks. Since 
time-dependent flows involve both the temporal and spatial errors, additional errors 
caused by mesh movement and dynamical refinement/coarsening are introduced 
into the numerical solution. The main source of error for moving-grid methods is 
the inability to satisfy the geometric conservation law in the presence of moving dis- 
continuities. One of the main problems associated with time-dependent local grid 
refinement methods is the need to interpolate. Note that any high-order interpo- 
lation across a strong discontinuity can cause a perceptible loss of spatial accuracy 
and can produce strong oscillations because of the Gibbs phenomenon. Another 
very important issue, particularly for unsteady problems, is a mutual influence of 
the temporal and spatial error components. 

8. CONCLUSIONS 

The accuracy of the adaptive grid redistribution and local grid refinement meth- 
ods is examined and analyzed for captured shocks. The grid refinement study using 
second- and fourth-order finite difference schemes based on the Lax-Friedrichs flux 
vector splitting is performed to solve the supersonic inviscid flow around a circular 
cylinder. The numerical calculations show that if the solution error is dominated 
by the first-order error component, the grid adaptation near the shock increases the 
error convergence rate towards the design order of the numerical algorithm used. At 
the same time, it has been shown theoretically and corroborated numerically that 
the design-order error component drastically increases because of the grid nonuni- 
formity. As a result, either clustering of grid points or local grid refinement near the 
shock improves the numerical solution accuracy only asymptotically on very fine 
meshes that are not reasonable for modern computers. Furthermore, it has been 
found that neither grid adaptation strategy practically reduces error in the pres- 
sure integral across the shock compared with that obtained on the corresponding 
uniform grid. From the present analysis it follows that grid adaptation is desirable 
if the first-order error component measured in the smooth part of the solution is 
much larger than the design-order error component. This inequality can be used 
as a grid adaptation criterion for captured discontinuities. 
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